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ABSTRACT 

We show that when time-reversible symplectic algorithms are used to solve periodic 
motions, the energy error after one period is generally two orders higher than that of the 
algorithm. By use of correctable algorithms, we show that the phase error can also be elim- 
inated two orders higher than that of the integrator. The use of fourth order forward time 
step integrators can result in sixth order accuracy for the phase error and eighth accuracy 
in the periodic energy. We study the 1-D harmonic oscillator and the 2-D Kepler problem 
in great details, and compare the effectiveness of some recent fourth order algorithms. 
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1 Introduction 



Symplectic integrators [1, 2, 3, 4, 5] preserve Poincare invariants when integrating classical 
trajectories. For periodic motion, their energy errors are bounded and periodic, in con- 
trast to non-symplectic Runge-Kutta type algorithms [6] whose energy error grows linearly 
with the number of periods[7, 8, 9]. Energy conservation alone suggests that symplectic 
algorithms are better long time integrator of classical motions. However, for periodic mo- 
tion, even symplectic algorithms are not immune from the linear growth of the phase error 
[7, 8, 9]. Whereas the energy error is the error of the action variable, the phase error is the 
error of the angle variable. Of the two, the phase error is even more important in determin- 
ing the long term accuracy of trajectories. For example, when symplectic algorithms are 
used to compute the Keplerian orbit, the elliptical orbit is easily seen to process. The pre- 
cession is of nearly constant radius. Since the semi-major axis of the ellipse is fixed by the 
initial energy, the constancy of the precession radius implies excellent energy conservation. 
Yet in spite of that, the precession itself implies that the trajectory is highly inaccurate. 
This orbital precession is a direct manifestation of phase error. Thus to preserve the long 
term accuracy of periodic trajectories, despite the primacy of energy conservation [10], one 
must seek to reduce the phase error directly. 

For periodic motion, the only error that matters is error that persists after one period [9]. 
A fundamental finding of this work is that, for periodic motion after one period, the energy 
error is at least (At)^ times that of the phase error, where At is the time step size used. 
Thus at small At the phase error is the dominant error governing the long term accuracy of 
periodic motion. Moreover, we show that the phase error of the symplectic corrector [11, 12, 
13, 14, 15] kernel algorithm is (At)^ times the phase error of other algorithms nominally of 
the same order. Recently, one of us[16] has made explicit the "correctability" requirement 
in deriving a correctable kernel algorithm. This criterion determines the optimal symplectic 
algorithms for solving periodic motion. The corrector algorithm has its origin in canonical 
perturbation theory[17]. It has been studied extensively [11, 12, 13, 14, 15] for its labor saving 
feature of only having to iterate the kernel algorithm. Here we draw the connection between 
symplectic corrector algorithms and the phase error in periodic motion. Much of our analysis 
is analytical rather than numerical, so that one can understand the result in a transparent 
way. We also found that forward time step symplectic algorithms[18, 19, 20, 21, 22] generally 
have much smaller phase errors than traditional algorithms with backward intermediate time 
steps[3, 5, 23, 24, 25]. 

In this work, we will analyze in detail the two fundamental prototypes of periodic motion: 
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the 1-D harmonic oscillator and the 2-D Kepler orbit. We are not interested in solving the 
harmonic oscillator per se, but only in using it as a vehicle for understanding the phase error 
and the working of our algorithms. It is only with such a simple model that we can show 
analytically how the phase error can be reduced by fine tuning the algorithm. To the extent 
that harmonic motion is the simplest periodic motion, this is clearly a necessary first step 
for proposing any scheme of phase error reduction. In the 2-D Kepler case, we demonstrate 
the usefulness of forward symplectic algorithms as compared to existing negative time step 
algorithms. For completeness, we begin with a brief review of the operator construction of 
symplectic algorithms, followed by a synopsis of symplectic corrector algorithms. In section 
5, we illustrate the basic idea of our analysis by showing how a second order algorithm can 
achieve fourth order accuracy in the phase error when solving the 1-D harmonic oscillator. 
In section 6, we repeat the same analysis for a class of fourth order forward algorithms. Error 
terms up to eighth order are computed by use of the Lie series [27] expansion. Beyond eighth 
order, the error terms can be determined by exactly solving the matrix model. All these 
are done analytically. We repeat the analysis for the Kepler problem in section 7. Here, 
we compare the phase error numerically for a number of recent fourth order symplectic 
algorithms. We summarize our conclusions in section 8. For the reader's convenience, some 
lengthy formulae and explicit calculations are given in the Appendix. 



2 Operator Factorization 

Symplectic algorithms can be derived most simply on the basis of operator factorization. 
(See the excellent review by Yoshida[2] and earlier references therein.) For any dynamical 
variable W{qi,pi), its time evolution is given by the Poisson bracket, and therefore by the 
corresponding Lie operator H associated with the Hamiltonian function H{qi,pi), i.e. 

dW_ ^ ^dWdH__dWdH_ 
dt ' ~ dqi dpi dpi dqi ' 

= mi.<^l.)w = HW. (2.2) 

^dpi dqi dqi dpi^ 

(Repeated indices imply summation). More generally, for any dynamical variable Q, we can 
define its associated Lie operator Q via the Poisson bracket 

QW = {W,Q}. (2.3) 

As we will see, this fundamental operator mapping underpins the entire development of 
symplectic integrators. 
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The operator equation (2.2) can be formally solved via 

W{t) = e^^W{0). (2.4) 

Symplectic algorithms are derived by approximating the evolution operator e* ^ for a short 
time in a product form. For Hamiltonian function of the standard separable form, 

//(q,p) =r(p) + y(q), with Tip) = ^pm, (2.5) 

the Hamiltonian operator (2.2) is also separable, 

H = f + V, (2.6) 

with first order differential operators T and V given by 

dT d _ d 
~ dpi dqi ^' dqi ' 

dqi dpi dpi 
Note that H, T and V individually satisfy the defining equality (2.3). 

The corresponding Lie transforms [27] e^'^ and e^^, are then displacement operators 
which shift qi and pi forward in time via 

q^q + £p and p^p + eF. (2.9) 

Thus, if e^^ can be factorized into products of Lie transforms e^-^ and e^^, then each 
factorization gives rise to an integrator for evolving the system forward in time. Most 
of the existing literature on symplectic algorithms is concerned with decomposing e^^ to 
arbitrarily higher order in the product form of 

N 

^e{f+V) ^J^^Uef^v.eV^ (2.10) 
i=l 

with a well chosen set of factorization coefficients {U^Vi}. In most cases, we will consider 
only the left-right symmetric factorization schemes such that either ii = and Vi = vjv-i+i, 
ti+i = tN-i+i, or Vat = and Vi = v^-i, U = tM-i+i- In either cases, the algorithm is 
exactly time-reversible, and the energy error terms can only be an even function of e. Such 

a symmetric factorizations is then at least second order. As first proved by Slieng[29], 
and Suzuki[30], beyond second order, decompositions of the form (2.10) must contain some 
negative coefficients ti and Vi. Goldman and Kaper[31] further proved that beyond second 
order, there must be at least be one pair of negative coefficients (tj, Vi). To circumvent this 
backward time step restriction[18, 19], one must faetorize the evolution operator in terms of 
operators T, V and the commutator [V^, [T, F]]. In this work, we will further demonstrate 
that these forward symplectic algorithms are also effective in reducing the phase error. 
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3 Symplectic Corrector Algorithms 

To see the relevance of symplectic corrector algorithms to periodic motion, we recapitulate 
some recent results[16]. Let Ta be a symmetric, approximate factorization of the short time 
evolution operator e^(-^+^\ 

N 

Ta = 1[ e*'^V*^^ = e'^^ , (3.1) 

i=l 

then the approximate Hamiltonian operator Ha must be even in e, i.e. 

H^ = f + V + e\eTTv[f, [f, V]] + evTv[V, [f, V]] ) + 0(£^) , (3.2) 

with error coefficients exTV^ ^VTV determined by factorization coefficients {ti,Vi}. Consider 
the similarity transformed propagator, 

= STaS-^ = Se'^^S-' = e^(^^^^'') = e'^'^ , (3.3) 

where the last equality defines the transformed Hamiltonian H'^. If now we take 

S = ex.p[eC] , (3.4) 

where C is the corrector, then the following fundamental result 

H'a = e'^HAe-'^ = Ha + e[C, Ha] + ^e^ [C, [C,Ha]] + -- - , (3.5) 

implies that 

H'^ = f + V + e\eTTv[f, [T, V]] + evTv[V, [f, V]] ) + e[C,f + ¥] + ■■■ . (3.6) 

One immediately sees that the choice 

C = ecTv[T,V], (3.7) 

would eliminate either second order error term with ctv = errv or ctv = evTV- More 
importantly, if Ha is constructed such that 

bttv = evTV , (3.8) 

then both error terms can be eliminated by the corrector. Thus for such an approximate Ta, 
the transformed propagator will be fourth order. This is the fundamental "correctability" 
requirement for correcting a second order T4 to fourth order [16]. In general, the corrector 
can be more complicated than the kernel algorithm Ta- However, when one iterates 7^, all 
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intermediate correctors cancel and only the initial and final corrector remains. For periodic 
motion, even the initial and the final corrector would have cancelled after exactly one 
period. Hence even if T4 is only second order, if it satisfies the correctability requirement 
(3.8), then its error after exactly one period would be fourth order! Thus among all second 
order algorithms, those that are "correctable", i.e. satisfy the the correctability requirement 
(3.8), would be two orders better. With a correctable algorithm, we will show later that 
the phase error is improved intrinsically even without applying the corrector. However, if 
the step size e is not commensurate with the period, one may step-over the minimum of 
the error function without knowing that it is there. In this case, it is essential to apply the 
corrector just prior to computing any observable. The advantage of a corrector algorithm 
is that for long-time integration, one usually only needs to apply the corrector sparingly at 
a few selected points in time. 

This correctability requirement can be generalized to higher order. At higher orders, 
Ha will have error terms of the form [T, Qi] and [V, Qi] where Qi are some higher order 
commutator generated by T and V. If Ha is of order 2n in e, then H'^^ can be of order 
2n + 2 if Ha^s error coefficients for [T,Qi] and [V^, Qi] are equal for each Qj. This is the 
fundamental corrector insight of [16]. In the following sections, we will demonstrate how 
this insight can be used to reduce the phase error in practical applications. 



4 The Modified Hamiltonian and Error Structure 

The distinct advantage of symplectic algorithms is not only that they preserve all Poincare 
invariants, but that their corresponding modified Hamiltonians and error structures can be 
systematically determined. This is of paramount importance when one seeks to understand 
the fundamental cause of an algorithm's error. To illustrate the approach, we begin by 
analyzing the simplest, first order factorization, 

where Ha is the approximate Hamiltonian operator 

Ha = H + ls[f, V] + ^e[f, [f, V]] - ^s[V, [f, V]] + .... (4.2) 

of the algorithm. This follows directly from Baker-Campbell-Hausdorff (BCH) formula. 
Thus the algorithm evolves the system according to the modified Hamiltonian Ha rather 
than the original Hamiltonian H. Nevertheless, the Hamiltonian structure of the system is 
preserved. As £ ^ 0, one recovers the original dynamics. Moreover, knowing Ha allows us 
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to determine the actual Hamiltonian function Ha which governs the algorithm's evolution. 
This can be done systematically by use of the Lie-Poisson bracket correspondence. To make 
this part of the discussion self-contained, we briefly summarize some pertinent results. 
From the fundamental defining equality (2.3), we can deduce Ha via 

HaW = {W,Ha}, (4.3) 

if we know how commutators of T and V transform back into functions under the operator 
mapping (2.3). By repeated applications of (2.3), we have 

[f, V]W = f{W, V} - V{W, T} , 

= {{W,V},T} - {{W,T},V} , 

= {W,{V,T}}, (4.4) 

where the last equality follows from the Jacobi identity 

{{W, V}, T} + {{T, W}, V} + {{V, T}, W}=0. 

Equality (4.4) implies the following correspondence between commutators of Lie operators 
and Poisson brackets of dynamical variables: 

[f,V]^{V,T} = -{T,V}. (4.5) 

There is thus a order reversal, or a simple sign change, in going from Lie commutators to 
Poisson brackets. (There is no such order reversal in the usual correspondence between 
quantum mechanical commutators and Poisson brackets.) This order reversal will only 
change the sign of odd-order brackets, as illustrated in the following examples: 

[V,[f,V]] {{V,T},V} = {V,{T,V}}, (4.6) 

[f,[V,[f,V]]] — > {{{V,T},V},T} = -{T,{V,{T,V}}}. 

Applying this to (4.3) gives, term by term, 

HaW = HW + ^e[f,V]W + ^e^[f,[f,V]]W-^e^[V,[f,V]]W + ..., 
{W,Ha} = {W,H} + {W,^e{V,T}} + {W,^e'{{V,T},T}}-... , (4.7) 

from which we can identify, 

HA = H-^s{T,V} + ^s''{T,{T,V}}-^e^{V,{T,V}} + ... . (4.8) 
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This general result merely transcribe expressions of Lie commutators into Poisson brackets. 
It is valid regardless of the form of the Hamiltonian. For the separable Hamiltonian (2.5), 
we have specific results 

dT dV 

{T.{T.n} = -j;j^=«V«. (4.10) 

Since T = T({pj}) and V = V^({^i}), there is no ambiguity about the meaning of subscripts 
on Tj or Vj. Also, since Tij = Sij, we therefore have, 

HA = H + ^epiVi + ^e^ PiVijp, + ViVi + . . . . (4.12) 

In general, the algorithm's approximate Hamiltonian is non-separable and more complicated 
than the original Hamiltonian. Similar expression has been given by Yoshida[2] in terms 
of Hp^, Hq.q., etc.. For a separable Hamiltonian of the form (2.5), one can certainly write 
Tj = Hp., and Vij = Hq.q., etc., but the latter is not more general than the former. If the 
Hamiltonian is not separable, Yoshida's expression suggests a degree of generality beyond 
that of the formalism. It is best to leave the form of the approximate Hamiltonian function 
in terms of Poisson brackets, which is then valid for all Hamiltonians. 

For higher order algorithms, the Hamiltonian operator corresponding to any left-right 
symmetric factorization is 

Ha = f + V + e'' [e^^y [f V] + e^^y [Vfv]) 

+ ( e^^^^y [f f3 V] + ey^^^y [V V] (4.13) 
+ errvTv [f{fvf] + ey,y,y [V{fvf]) + ... , (4.14) 

are coefficients specific to a particular algorithm and where we 
have used the condensed commutator notation [T^F] = [T, [T, V\\. Note that for symmetric 
decompositions, one has only even order commutators and the Lie-Poisson correspondence 
is trivial. In terms of similarly condensed Poisson brackets, {T'^V} = {T, {T,V}}, the 
Hamiltonian function can be read off by inspection. 

Ha = T + V + e'' [e^^y{T''V} + ey^y{VTV}) 

+ ( e^^^^y {T V} + Cy^^^y {V V} 

+ eTTVTv{T{TVf] + ey,y,y{V{TVf]) + ... . (4.15) 
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For the separable Hamiltonian (2.5), these higher brackets are: 

{TT^V} = PiPjPkPiVijkl , 

{VT^V} = -SpiPjVijkVk , 

{TiTVf} = -2pi{VikjVk + VikVkj)pj, 

{V{TV)'} = 2ViVijV,. (4.16) 

The results in this section will allow us to analyze any symplectic algorithm from second to 
sixth order. Beyond sixth order, the number of Lie and Poisson brackets proliferates and 
other means of determining the Hamiltonian error terms may be more efficient. 



5 Harmonic Oscillator: Second Order Integrator 

To illustrate some of our key ideas in the simplest context, we will begin our study of the 
phase error with the second order factorization scheme 

T2{s,a) = e^^^e^^^e^^^, (5.1) 

with Vi given by 

Vi = V + ae^ [V, [f, V]] . (5.2) 
Classically, this Lie commutator produces a modified force [19] 

[V, [T,V]] = 2F,-1— = V,|F|2_ , (5.3) 
dqj dpi dpi 

resulting in the following more general second order symplectic integrator 

1 

qi = qo + 2 £ Po > 



Pi = Po + £ 



F(qi)+a£^V|F(qi)|^ , (5.4) 



1 

q2 = qi + 2 ^ Pi • 

Here, (qo,Po) a-^d (q2)Pi) are the initial and final states of the algorithm respectively. 
The introduction of the gradient term with parameter a will allow us to satisfy the cor- 
rectability criterion in its simplest setting. When applied to the 1-D harmonic oscillator 
with Hamiltonian 



2 ' 2 

the force gradient is just 



Hiq,p) = ^ + l-u;'q\ (5.5) 



F{q) = -u^q V,\F{q)\^ = 2u;^q. (5.6) 
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^TTVTV AOn Oyl ^ ' ^VTVTV ion ' (^■^) 



For the standard Hamiltonian, the approximation Hamiltonian operator for any sym- 
metric factorization is given by (4.14). The non- vanishing error coefficients corresponding 
to algorithm (5.1) are just 

1 _ 1 

24' ^^^^-""T^ 

11 _ 1 1 

480 ~ 24" ' ""^^^^^ ~ 120 ~ 6' 
The Hamiltonian function is then as given by (4.15). For the harmonic oscillator as defined 
by (5.5), we have Vij = u^Sij, Vijk = 0, {TT^V} = 0, {VT^V} = and non- vanishing 
brackets, 

{T,{T,V}} = uV, 
{V,{T,V}} = -a;V, 
{T{TVf} = -2a; V, 

{V{TVf} = iJ^i. (5.9) 

Notice the clear separation between the contributions of the algorithm, which are the error 
coefficients, and that of the physical system, which are the Poisson brackets. The final form 
of the Hamiltonian function due to algorithm (5.4) is therefore, 

-2uj^e^ {^e^^y^yP^ - ey^y^yUp q^^ + ... , (5.10) 

= 7r—P^ + \k*q^. (5.11) 
2m* ^ 2 ^ ^ ' 

Thus the oscillator being evolved by the algorithm is one with an effective mass and spring 
constant, 

m* = m*{e) = {I + 2 e'^ ixi"^ e^^y - 4.e'^ uj^ e^^y^y + . . .)-^ , (5.12) 
k*=k*{e) = {l-2e'^ up ey^y +Ae'^uj^ey^y^y + ...)J^ , (5.13) 

from which one can deduce the approximate angular frequency 

UAie) = J^. (5.14) 
V m* 

The phase error is simply related to the fractional deviation of the the approximate angular 
frequency from the exact frequency: 

A(j)={uA-co)T = 2Tr{—-l). (5.15) 
10 



This is the fundamental thrust of our analysis: tracking the phase error of the algorithm 
back to its factorization coefficients. Observe now that from (5.12) and (5.13), we have 



Me) = a;^(l + 2£2a;2e^^^+ ...)(l-2£2a;2e^^^+ ...), (5.16) 



CO 



l + e^u\e^^^-e^^^) + 0{s^) . (5.17) 



In general, the approximate frequency is second order in error, as befitting a second order 
algorithm. However, if the correctability criterion e^^,^ = e^^,^ is satisfied, then uja is 
fourth order. Moreover, if the algorithm is originally fourth order with 6^,^,^ = eyj,y = 
then satisfying 6^^,^^,^, = e^j,^,^,^ would make loa sixth order. Thus an nth algorithm can 
have an (n+2)th order phase error if its error coefficient satisfies the correctability criterion. 
This is the key connection linking the phase error with correctable algorithms. (Note that 
by making e^^,^ = Sy^y (but not zero) and 6^^^^^^ = Sy^y^y, would not make the phase 
error sixth order.) 

With only one free parameter presently available, we can only set e^^y = eyj.y = — ^ 
with the choice 

« = ^ , (5.18) 

thus making loa fourth order. This particular value corresponds to the well known propaga- 
tor first derived by Takahashi and Imada[26] for computing the quantum statistical trace [26] 
to fourth order. The same factorization scheme, interpreted as symplectic corrector algo- 
rithm (5.4), has also been used by Lopez-Marcos et al.[13, 14] and Wisdom et al.[ll] for 
solving classical and celestial dynamical problems. With this choice of a, the coefficient of 
the fourth order frequency error is, from (5.12), (5.13) and (5.8), 

,(4) 



= lim 



1 f'^A 



£4 V CO 



2iO [eyj:yj:y " ^ rj^rj^y " ^TTVTV) ~ ^ ^JQQ ' (5.19) 



To gauge the relative importance of this phase error, let's compare it to the energy error 
after one period. Since it is the modified, or approximate Hamiltonian that is conserved by 
the algorithm, i.e. 

HA{q,p)=HA{qo,Po), (5.20) 
the energy after one period T = 27r/a; can be expressed as 

Hiq^,p^) = H{qo,po)+e^AH^^\e^)+e^AHi^\e^) + e'^AHi^\£^) + 0{£'') . (5.21) 
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Prom (5.10), we have in particular, 

Ai??^(£') = - u;-" (e^^^ip-'-pD- e^,^u\q^-ql))\^^^, (5.22) 

In order to compute these energy deviation errors, we must solve for p{t) and q{t) according 
to Hamiltonian Ha'- 

/ q{t; e)\ / cos{uJAt) (m*^^)"-^ sin(a;Ai) W ^ ^-^ 24) 

Since m* and loa are e^-dependent, each function AH^^^e"^) contains further dependence 
on We now define the constant energy error coefficients Ej^"^ via 

H{q^,p^) - H{qo,po) = AEt = E%^ + E%^ + E%^ + 0{e^) , (5.25) 

where for example, we have 

= A<)(0) + Ai7(,')'(0), 
eP = AHP{0) + AHi^^'iO) + ^AH^^^"{0), 

= AhP{0) + AhP'{0) + ^AH^^'^"{0) + ^AH^^^"'{0). (5.26) 

Here, the prime denotes derivative with respect to . Prom the form of each AhS^ (e^^., 
since £ = implies that uja = oj, p(T) = po and g(T) = qo, we must have 

Ai?(T^(0) = 0, (5.27) 

and therefore 

E%^ = . (5.28) 

Thus for periodic motion, despite the fact the algorithm is only second order, the energy 
error is actually fourth order after one period. 
The fourth order energy error is given by 



= Ai7?)'(0) = -2a;2 (e^^^pTPx' - e^^^a^V 9t)| > 

= Attu^po qo {e^^y - ey^y){erj.^y + Cy^y) , (5.29) 



where we have used 



g'(T; 0) = - w^(0)T and p'{T- 0) = -oo qo u;'a{0)T , (5.30) 
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and from (5.17), 

'^a(0)T = 2-Ku?{e^^^, - Cy^y ) . (5.31) 

The fourth order error now vanishes if the algorithm satisfies the correctabihty criterion 
ej.j.y = eyj.y . TYuis for a correctable second order algorithms, after each period, the phase 
error is fourth order and the energy error is sixth order. 

Since the factor (5.31) is common to all first derivatives (in e^), we conclude that for 

^TTV ^VTV 

AhP'{0) = 0. (5.32) 

Hence for e^^,^ = Cy^y , the sixth order energy error can be now computed as 

E^^ = ^AiJ^^"(0), (5.33) 



2 



VTV ) 



4:Trp() qoU) {^ttv ~^ ^vtv) "^{^ttvtv ^vtvtv) ~^ ^ttv ^vtv 

Po qo ■ (5.34) 
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2160' 

The above calculation demonstrates the general property of the energy deviation error 
after one period. For correctable algorithms, the first two terms in the error expansion (5.26) 
vanish identically, which means that to compute e!^'\ one need not know the explicit form 
AhI^\£'^). However, in order to compute Aii'^-*"(0), one must know m*(£^) and loa{s'^) 
accurately to 0{e^), which means knowing the fourth order Hamiltonian error function, 
or Ai7^^(£^). Thus although (5.33) makes no reference to Ai?^''(e^), one must know it 
implicitly. Similarly, eI^^ can be computed from AH^rp^e^) and AH^\e'^) via 

4«) = ^Ai?¥)"(0) + ^Ai?¥)'"(0) . (5.35) 

However, in order to compute AiJ^^'"(0) one must know AH^'^\e'^) correctly to O(e^). This 
would again require knowing the sixth order error Hamiltonian or AH^^\e^). In general, 
e!^^ can be compute two orders beyond the accuracy of knowing the Hamiltonian. 

To summarize, for a second order algorithm, the energy after one period is automatically 
fourth order in e (= At). If the algorithm is correctable, then the energy error is sixth order. 
For special initial conditions po = or go = 0, by solving the algorithm exactly in the case 
of the harmonic oscillator [28], one can show that the energy error is actually tenth order. 
This last error reduction only occurs for the harmonic oscillator. Nevertheless this further 
emphasizes that the energy error after one period is not a very good gauge of any integrator's 
accuracy. On the other hand, the phase error, as reflected in the fractional change of the 
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oscillator's angular frequency, can at most be fourth order and is a much more stringent 
and discriminating benchmark. 



6 Harmonic Oscillator: Fourth Order Forward Integrators 

Beyond second order, all symplectic algorithms of the form (2.10) must have some neg- 
ative intermediate time steps[29, 30, 31]. This means that at some intermediate time, 
the algorithm is moving the phase trajectory backward in time. For classical mechan- 
ics, which is timc-rcvcrsiblc, these negative time steps arc harmless. However for solving 
timc-irrcvcrsiblc problems, such as the diffusion or Fokkcr-Planck equation, backward time 
step evolution is not possible. These systems can only be solved by forward decomposition 
algorithms, with all positive, even intermediary, time steps. Some fourth order forward 
algorithms have been derived recently for solving a variety of time- irreversible [32, 33], and 
time-reversible[19, 21, 22] equations, both with excellent results. Beyond second order, 
purely forward time steps are possible only if one include the commutator [V, [T, T]] in 
addition to operators T and V in the factorization process. In this work we will apply these 
fourth order forward algorithms to study the phase problem of periodic motion. In this 
section, we further generalize our study of the harmonic oscillator by use of these fourth 
order forward algorithms. 

Chin and Chen[21, 22] have introduced a family of fourth order forward algorithms 
4ACB parametrized by a parameter to. We use here a slightly generalized form by mul- 
tiplying the central commutator by 1 — a and adding a/2 times the commutator to each 
potential operator on each side. The resulting algorithm has the operator form 



gto £ T^Vl £ Vl 1 £ Tgj;2 £ V2 1 £ T^Vl £ Vl gto £ T 



(6.1) 



where 



Z Vl 

t> + (l-a)-£' [V,[f,V]], 



12 L 1 - 2to 6(1 - 2to)^ 



(6.2) 



(6.3) 



and 




(6.4) 



The corresponding forward symplectic integrator can be read off directly as 



qi = Qo + stopo, 
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t'iF(qi) + -^Xoe'V|F(qi] 



Pi = Po + e 
q2 = qi + £ Pi , 

V2 F(q2) + (1 - a) no e^V\F{q,2)f\ , (6.5) 



P2 = Pi + £ 

qs = q2 + £tiP2, 

P3 = P2 + £ 

q4 = q3 + eioP3, 



^iF(q3) + -woe V|F(q3)r 



where (qo , po) and (q4 , P3) are the initial and final states of the algorithm respectively. The 
parameter a can be changed from to 1, but there is really no restriction on its range. When 
applied to the harmonic oscillator, the parameter a can be used to correct the algorithm to 
sixth order. The parameter to can be varied from to ic = ^(1 — ~ 0.21. For to = 0, 
the final force evaluation can be reused at the next iteration, thus eliminating one force 
evaluation. At the upper limit of to = tc, V2 = 0, also eliminates one force evaluation. For 
to > tc, V2 becomes negative, and the algorithm ceases to be a forward algorithm. 

Our analysis of the second order algorithm can now be repeated verbatim for the fourth 
order case. The approximate Hamiltonian operator corresponding to any symmetric fourth 
order algorithm is of the form, 

HA = f + V + e^ ( e^^^^^[ff^V]+e^^^^y[Vf^V] (6.6) 

+e,,^,^ [f{fvf] + e^,^,^ [V{fvf] ) + 0(£6) . 

For the harmonic oscillator, [T'^y] = 0, and the first two error term vanishes identically. 
The evaluation of the last two error coefficients for the family of fourth order algorithm 
(6.5) is non-trivial and is given Appendix A. The corresponding Hamiltonian function, 
after recalling the Poisson form (4.15) and brackets (5.9), is 

HA{q,p) = Y + {^TTVTVP^ - (^VTVTV^'^ (f) + ■■■ ^ (6-7) 

^ P^ + W<1^^ (6-8) 



2m* 2 
with 

m* = m*{e) = {1 - Ae"^ e^^^^^ + . . .)-^ , (6.9) 

k* = k*{e) = u;^{l + ie^u;'^ey^^^^+ ...), (6.10) 

and approximate frequency 



ua{£) = CO J {l + Ae^io^e^^y^y +..:){! -Ae^u'^e^^^^y + ...), (6.11 



to 



1 + 2e^ io\ey^^^^ - e^^^^y) + Oie"^) . (6.12) 
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Again, one immediately sees that if the sixth order correctabihty criterion 



is satisfied, then ua will be sixth order. Note that now we have 



(6.13) 



u'aT 



= 0, 

£=0 



^'a^ = ^T^i^ {errvTV - (^VTVTv) ^ (6-14) 

where primes still denote derivative with respect to e^. The conservation of HA{q,p) again 
implies that the energy deviation after one period can be expressed as 

H{q^,p^) = H{qo,po) + £^ AH^^\e'') + AH^ie'') + AH^^\e'') + 0(6^") , (6.15) 

with 

AHi^\E'') = (e,,^,^ {p'-pD - e^,^,^ a;^ (g^ - q^)) | . (6.16) 
The constant energy error coefficients Ej!'^ defined by 

H{q^,p^) - H{qo,po) = AEt = E^^^ + E%^ + E^ + e^^ E^^""^ + 0{e^^) , (6.17) 

are now of the form 

4") = AH^4\(}), 

eP = AiJ?)(0) + A//^')'(0), 

eP = Ai?(,^)(0) + Ai7f (0) + iAi7^)"(0), 

= AH^^''\0) + AHi^^'{0) + ^AHP"{0) + ^AH^^^"'{0). (6.18) 
Now, because of (6.14), for e^^^y^ = eyj,yj,y, not only we do have AH^^\o) = 0, but also 

AhP'{0) = and AH^^'^" (0) = . (6.19) 

This implies that 

=4^) =4^) =0, (6.20) 
and the first non- vanishing energy error is tenth order, 

4° = iAif^)'"(0). (6.21) 

However, as noted in the last section, in order to compute this, one must determine the 
sixth order error Hamiltonian. 
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Due the complexicity of the algorithm, these higher error terms are difficult to compute 

by Lie series. However, they can always be computed using the matrix method[28]. For 
brevity, we will skip over the details and just report the final results. 

We have shown earlier that the fourth order phase error term will vanish if ej,j,yj.^ = 
^vTVTv -^'^'^ ^ given value of to, this criterion can now be satisfied by a specific choice of a 
given by a = aito) in (A. 12). Using this functional form to eliminate a in terms of Iq, the 
sixth order error term = /(^o) scaled such that a; = 1, is plotted in Fig.l. 



1.2x10' 




-4.0x10 



0.00 



0.05 0.10 



0.15 



0.20 0.25 



Figure 1: The sixth order angular frequency error as a function of the algorithm's parameter 

h. 

Within the forward range of < < 0.21, the sixth order frequency error has a 
minimum of value 

^(6) 

= 7.718621317057857 x 10-^a;^ (6.22) 



at to = 0.12129085056575276, and a pole at = 0.13882413776781183. Note that outside 
of the forward range, the error can actually vanish at to = 0.24265927253055103. 
The eighth order energy deviation error after one period is 

A^;?^ = 16 TT ^^^^ _ ^2 ^^^^ ) ^ (g_23) 

which 

second order case. 
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Thus for a corrected fourth order algorithm, the first non-zero energy deviation error is 
tenth order. This is plotted in Fig.2 scaled such that uj = go = po = 1. 




Figure 2: The tenth order energy deviation error after one period as a function of the 
algorithm's parameter to- 

Within the forward range of < to < 0-21, the tenth order energy deviation error has a 
minimum of value 



= -1.3398713813012635 x W'^ uj^^ qqPq , (6.24) 

mm 



at to = 0.12482248354859667, and a pole at to = 0.13882413776781183 (same as in the 
frequency case). In both cases the error term vanishes at the same value, i.e. to = 
0.24265927253055103, outside of the forward range. (Note also that this error term vanishes 
for special starting value of po = or go = 0. It can be shown that for either po = or 
qo = 0, the first non-vanishing energy error term is 16*'* order, again demonstrating that 
the phase error dominates overwhelmingly over the energy error.) 

7 The 2-D Kepler Problem 

In light of our previous discussion, for long term trajectory simulation, one must judge all 
symplectic algorithms on how well they minimize the phase errors rather than the energy 
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error. In this section, we will examine Keplerian motions in 2-D defined by the Hamiltonian 

Here, our analysis of fourth order algorithms will not be as extensive as in the harmonic 
oscillator case because the approximate Hamiltonian 

HA = f + V + e^ ( e^^^^^ [ff^V] + e^^^^^ [Vf^V] (7.2) 

can no longer be solved analytically. The operator [T^V] ^ and while we can still 
force ej,j,yj,v = &vtvtv ^ ™ harmonic oscillator case, we have no way of ensuring 
that ej,j,j,j,y = eyj,j,j,y. Currently, there are no known fourth order forward symplectic 
algorithms that can be corrected to sixth order. Nevertheless, identical analysis as in the 
harmonic oscillator case shows that 

E^^ = AH^}\0) = , (7.3) 

and the energy error after one period must be at least sixth order. Thus if fourth order 
algorithms are used to solve Keplerian orbits, it is more fitting to examine their fourth order 
phase errors instead. 

For two-dimensional motion, there are two basic phase angles associated with the two 
sets of canonical variables {qi,pi) and {q2,P2)- A convenient measure of these phase errors 
is the precession error of the orbit in the {qi,q2) plane, which can be tracked[20] by the 
rotation of the Laplace-Runge-Lenz (LRL) vector 

A = pxL-q. (7.4) 

In the above definition, L = q x p, is the angular momentum vector. 

To sec how various algorithms compare, we first plot the fourth order energy error 
function defined by 

H,{ci{t)Mt)) = lim -^[£;(q(t),p(t)) - Eo] , (7.5) 

in Fig.3. 
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Figure 3: The energy error at half a period for an eccentricity of 0.9 



Note that this is an intrinsic function characteristic of each algorithm independent of the 
step size. We compute this function by finding the energy deviation from the initial energy 
along the orbit and then dividing it by e^. As £ gets smaller and smaller, this function 
converges to its limiting form. The functional form is basically unchanged for e < T/3000, 
where T is the period of the Keplerian orbit. All results shown in Fig. 3 are computed with 
e = T/5000. 

Since wc have shown that -E'(q(T), p(T)) — Eq = 0{e^), H4 vanishes exactly after one 
period. Thus each of energy error curve of Fig. 3 reverts back to zero at t = T. This is 
a characteristic behavior of all symplectic algorithms. Non-symplectic Rungc-Kutta algo- 
rithms do not have this property and their energy deviation error accumulates rather than 
vanishing after each period. However, even for symplectic algorithms, the energy deviation 
error is non-vanishing at other times. Here, due to the high eccentricity (e = 0.9) of the 
orbit, the energy error is at a maximum near mid-period. Algorithm Chin-C (C), is the 
forward algorithm (6.1) with to = 1/6 and a = 0, first derived in [19]; Blanes-Moan (BM) 
is an algorithm recommended in McLachlan and Quispel's review[5]; Omelyan et a/.[25](0) 
is a recent alternative forward algorithm that uses the same force gradient defined by (5.6); 
McLachlan [3] (M) is a greatly improved version of the first fourth order Ruth- Forest [23] 



20 



algorithm. With the exception of M, all algorithms have comparable error height at mid- 
period. Note however that BM requires six force evaluations, M uses four force evaluations, 
O uses four force plus four force- gradient evaluations, but C uses only three force and one 
force-gradient evaluation. Algorithm M's error height reaches up to 14, which is more than 
twenty times higher. This is rather surprising, since algorithm M works very well in solving 
quantum mechanical [21, 34] and three-body [22] problems. 

In Fig. 4, we track the rotation of the LRL vector during orbital motion. 




Figure 4: The precession deviation error after half a period for eccentricity 0.9 with starting 
point q = (10, 0) and p = (0, 0.1) 



If the orbit is exact, the LRL vector is a constant vector pointing along the semi-major 
axis of the orbit. If the orbit precesses, then the LRL vector rotates accordingly. At any 
point in the orbit, the angle of the LRL vector is given by 

-Ay{ty 



e{t) = tan 



e%(t) + e''06(t) + 



and from which one can extract the fourth order angle error function via 



(7.6) 



(7.7) 



Again, this intrinsic function is computed in the limit of small e. We have checked that 
it has indeed converged to its limiting form for e = T/5000. Since the orbit precesses 
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the most when the particle is closest to the attractor, the LRL vector rotates measurably 
only during mid-period. It is constant before, and remained constant after the mid-period. 
Thus the rotation after one period is essentially the same as the rotation shortly after mid- 
period. Note that this (phase) angle error do not revert back to zero after each period, but 
accumulate after each period even for symplectic algorithms regardless of order. Thus the 
only way to minimize this phase error is to make it as small as possible. From Fig. 4, we see 
that algorithm C's rotation angle after mid-period in nearly an order of magnitude smaller 
that that of cither BM or O. The actual values after one period are: 0.0076, -0.0692, -0.1466 
respectively. Algorithm M's rotation function reaches down to —2.5, which is an order of 
magnitude greater than that of BM and O and two orders of magnitude greater than that 
of C. We did not bother to plot it. 

Since parameters to and a are at our disposal, we can further optimize the family of 
algorithm (6.1) to reduce the rotation error. The resulting optimal choice is shown in Fig.5, 
with to = 0.166160 and a = 0. The angle error after one period is further reduced by a 
factor of five from 0.0360 to 0.0077. 




Figure 5: The precession deviation error after half a period for eccentricity 0.936 with 
starting point q = (10, 0) and p = (0, 0.08) 

While one can optimize the family of algorithm (6.1) for any one specific problem, or at 
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one eccentricity, it is of greater value to devise an optimal algorithm for solving a general 
class of problems. For the Kepler problem, all possible shapes of closed orbits are spanned 
by the eccentricity; it is thus more desirable if one can devise an optimal algorithm for all 
values of the eccentricity. In Fig. 6, we plot the LRL rotation angle after one period as a 
function of the orbit's eccentricity, as determined by different initial conditions. 
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Figure 6: The precession deviation error for highly eccentric orbits. 

Most algorithms work well for orbits of low eccentricity and the rotation angle is cor- 
respondingly small. We therefore compare algorithm at e > 0.9 . At e = 0.95, the angle 
error values for M, BM, O and C are respectively -166.1870, -4.8865, -10.4470 , and 0.1244. 
Algorithm C's angle error is orders of magnitude smaller than other algorithms. 

In Fig.7, we again show that a better algorithm can be devised from the family of 
algorithms (6.1). The choice of a = (only one force-gradient), and to = 0.166160 (only 
slightly below the canonical value oito = 1/6), produces an algorithm with uniformly small 
phase error up to e = 0.95 . At e = 0.95 the angle error value for Opt-C is -0.00357, 
compares to C's value of 0.12363. 
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Figure 7: The precession deviation error for highly eccentric orbits. 

8 Conclusion 

In this work we showed that for periodic motion, the energy error after one period is 
generally two orders higher than that of the algorithm. If the algorithm is correctable, the 
phase error can also be reduced two orders higher. The use of fourth order forward time 
step integrators can result in sixth order accuracy for the phase error and eighth accuracy 
in the periodic energy. By generalizing the recently discovered one-parameter family of 
fourth order symplectic algorithms [21], we can minimize the energy and phase error to even 
higher order. The results of this study provides a direct verification of Chin's correctability 
criterion [16] for correcting a symplectic algorithm to higher order. In particular, wc showed 
that the correctability criterion is superior to the conventional wisdom of minimization of 
the sum of squares of error coefficients. The most important conclusion of this work is 
that for periodic motion, the phase error is a more discriminating gauge of an algorithm's 
effectiveness than the energy error. 

As a more important application of the phase error analysis, we track the orbital pre- 
cession angle of the 2D Kepler problem by monitoring the rotation angle of the Laplace- 
Runge-Lenz vector [20]. By comparing with various recent fourth order algorithms, we 
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demonstrated the uniqueness of forward symplectic algorithm in minimizing the phase er- 
ror of this important class of celestial mechanics problems. 
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Appendices 

A Fourth Order Error Coefficients 

The error coefficients of the fourth order forward algorithm (6.1) can be computed in terms 
of algorithm's factorization coefficients via a Mathematica program[32]. They are: 



= 2{to + ti), 

Cy = {2VI+V2), 



1 



tj (-4 vi + V2) + tl{2vi+V2) + 2to h (2 vi + V2) 



6uo-to (2 vi + ^2)^ + h {2vl + 2 vi V2 - vl) 



1 



Itl {to + Ui){2vi+V2) 



360 

+ t\ (4 to + h) {7v2 -16vi) + 6 t^tj {4:vi + 7v2] 
1 

90 



24 {to + 3h){2vi + V2y 



-Qtotj {6vl +V1V2- vl) + 4 i8vl -7viV2 + 2i;|) 



1 

60 



4{2vi + V2f + tl (lO{3a-l)uo + h{-lQvl + AviV2 + vjfj 



(A.l) 
(A.2) 
(A.3) 

(A.4) 
(A.5) 
(A.6) 

) 

(A.7) 



+ tl (-lOtio + ii i2vl + 2viV2 + Svl)) 
+ to tl (-20 no + tl (12 vf + 2viV2 + 3 V2)) 

[2tl{2vi+V2f -4toi2vi+V2) (5uo + ti{vl + viV2-vl)) (A.8) 
+ ti (l0no(2?;i + (3a-2)t;2)- ti {Avf + vlv2 + Sviv^ - 2v^)'^ 



In order for the algorithm to be fourth order, wc must have Cj, = Cy 
^vTv ~ 0- These four constraints can be satisfied by 

1 1 

tl=t2 = --to, ts=to, ^l=^3=g^^_2^^p, 



1 and e„ 



(A.9) 
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1 

= 1 - (Ui + V3) , Uq = — 



1 1 



(A.IO) 



1-2*0 6(l-2to)^ 

This is the family of fourth order algorithms (6.1) with parameters to and a. For the 
harmonic oscillator, 6^^^^,^ and e^^y^^ vanish identically. A fourth order algorithm can 
be corrected to sixth order if one can set e^^^y^ = e^^^^^. Substituting (A. 9) and (A.IO) 
into (A. 7) and (A. 8), gives e^^^r^y and eyj.yj.y as functions of the parameters to and a, i.e. 

1 + 5 a - 12 to (1 + 5 a + 20 a to (-1 + to)) 



(A.ll) 



^^^^^ 2880(1 - 2 to) 

_ l + 10a-6to(3 + 30a-to(9 + 210a + 8to(l-85a-3to(l-40Q; + 20Q;to)))) 
''^'^^''^ ~ 4320(1 - 2 to)4 

Solving for e^^yj,y = ey^yj,y determines a as a function of to: 

^ ^ l + 6<o(-3+ i/o(6 + /o(-23 + 2i/„))) 

5(l-12to(l-2to)2)(l-6to(l + 2to-4t2)) ■ ^ ' > 

However, there exists no real solution of the parameters for which both, e^^^j,^ and eyj.y^y 
can be set to zero, i.e. , we can have an algorithm that is correctable to sixth order, but 
not a real sixth order algorithm. 
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